1. Wykorzystanie analizy regresji do określenia zależności między ogólną liczbą ludności a punktami adresowymi.

maxpop = max(c(out_model_punkty_adresowe$TOT, out_model_punkty_adresowe$EST_POP))

p1 <- ggplot(data = out_model_punkty_adresowe) +
  geom_sf(aes(fill = TOT)) +
  scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop)) + 
  labs(title = "Ryc.1 Mapa oryginalnych wartości liczby ludności") + 
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

p2 <- ggplot(data = out_model_punkty_adresowe) +
  geom_sf(aes(fill = EST_POP)) +
  scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop)) + 
  labs(title = "Ryc.2 Mapa estymowanych wartości liczby ludności") + 
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

print(p1)

print(p2)

p3 <- ggplot(data = out_model_punkty_adresowe) +
  geom_sf(aes(fill = RES)) +
  scale_fill_gradient2(name = "Błąd estymacji", low = "blue", high = "red", limits = c(min(out_model_punkty_adresowe$RES), max(out_model_punkty_adresowe$RES))) + 
  labs(title = "Ryc.3 Mapa różnic między wartością obserwowaną liczby ludności a jej estymacją") + 
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

print(p3)

out_model_punkty_adresowe$RES_reclass = case_when(
    out_model_punkty_adresowe$RES < 0 ~ "bledy ujemne",
    out_model_punkty_adresowe$RES == 0 ~ "brak bledu",
    out_model_punkty_adresowe$RES > 0 ~ "bledy dodatnie",
  )

p4 <- ggplot(data = out_model_punkty_adresowe) +
  geom_sf(aes(fill = RES_reclass)) +
  scale_fill_manual(name = "Typ błędu", 
                    values = c("bledy dodatnie" = "red", "brak bledu" = "white", "bledy ujemne" = "blue")) +
  labs(title = "Ryc.4 Mapa zreklasyfikowanych błędów modelu") +
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

print(p4)

meanerror = mean(out_model_punkty_adresowe$RES)
rmse = sqrt(mean((out_model_punkty_adresowe$RES)^2))

print(paste("średni błąd estymacji modelu wynosi", round(meanerror,2), "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))
## [1] "średni błąd estymacji modelu wynosi 3.67 natomiast pierwiastek średniego błędu kwadratowego wynosi 30.44"

1.1. Ocena modelu

Maciejo,
YOU KNOW WHAT TO DO

2. Wykorzystanie analizy regresji do określenia zależności między ogólną liczbą ludności a liczbą mieszkań.

2.1 Liniowy model zależności między liczbą ludności a liczbą mieszkań w siatce 1km

# wykres rozrzutu
ggplot(pop1km, aes(x=TOT, y=N_BUDYNKI)) + geom_point() +
  coord_fixed(ratio = 1) + 
  xlab("Liczba ludności") + 
  ylab("Liczba mieszkań") + 
  geom_smooth() + 
  geom_abline(color = 'grey') + 
  theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

# model liniowy
model_liniowy <- lm(TOT ~ N_BUDYNKI, data = pop1km)
plot(model_liniowy)

W jakim stopniu liczba ludności (TOT) jest wyjaśniana przez liczbę mieszkań?

2.2. Jakie są statystyki reszt? (residuals)

resid_interact(model_liniowy)
## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## ℹ The deprecated feature was likely used in the ggplot2 package.
##   Please report the issue at <https://github.com/tidyverse/ggplot2/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
summary(model_liniowy$residuals)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -98.6520  -2.5883  -0.7094   0.0000  -0.5883 174.7735
mean(model_liniowy$residuals)
## [1] 3.738246e-16

3. Ocena dopasowania modelu (wykresy diagnostyczne, identyfikacja wartości odstających).

Wartości odstające można zidentyfikować wykorzystując wykresy diagnostyczne, a w szczególności wykres Residuals vs. Leverage.

autoplot(model_liniowy) + theme_classic()

#dopasowanie modelu
smr_model_lm <-summary(model_liniowy)

c( R2 = smr_model_lm$r.squared, R2_adj = smr_model_lm$adj.r.squared)
##        R2    R2_adj 
## 0.7019419 0.6990481
#~70% zmienności zmiennej TOT jest wyjaśniona przez zmienność zmiennej N_BUDYNKI.

# Określenie dowolnego przedziału ufności dla współczynników modelu
confint(model_liniowy, level = 0.99)
##                 0.5 %   99.5 %
## (Intercept) -7.022728 8.441527
## N_BUDYNKI    3.275627 4.603257
outlierTest(model_liniowy)
##     rstudent unadjusted p-value Bonferroni p
## 21  9.742598         3.0248e-16   3.1760e-14
## 25  6.813338         6.7995e-10   7.1395e-08
## 34 -5.176442         1.1389e-06   1.1958e-04
ggplot(pop1km, aes(x=TOT, y=N_BUDYNKI)) + 
  geom_point() +
  geom_point(data = pop1km[c(21, 25, 34),], aes(x=TOT, y=N_BUDYNKI), color = "red", size = 4) +
  coord_fixed(ratio = 1) + 
  xlab("Liczba ludności") + 
  ylab("Liczba mieszkań") + 
  geom_smooth() + 
  geom_abline(color = 'grey') + 
  theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

par(mfrow = c(2, 2))
plot(model_liniowy)

meanerror = mean(model_liniowy$residuals)
rmse = sqrt(mean((model_liniowy$residuals)^2))

print(paste("średni błąd estymacji modelu wynosi", meanerror, "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))
## [1] "średni błąd estymacji modelu wynosi 3.73824648485021e-16 natomiast pierwiastek średniego błędu kwadratowego wynosi 26.25"

PYTANIA Czego dowiadujemy się z wykresów diagnostycznych? Czy model jest dobrze dopasowany? Czy spełnione są założenia regresji liniowej?

4. Model po usunięciu wartości odstających

4.1 Model regresji

# Usunięcie wierszy 21, 25 i 34
pop1km_n <- pop1km[-c(21, 25, 34), ]

# wykres rozrzutu
ggplot(pop1km_n, aes(x=TOT, y=N_BUDYNKI)) + geom_point() +
  coord_fixed(ratio = 1) + 
  xlab("Liczba ludności") + 
  ylab("Liczba mieszkań") + 
  geom_smooth() + 
  geom_abline(color = 'grey') + 
  theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

# model liniowy
model_liniowy_n <- lm(TOT ~ N_BUDYNKI, data = pop1km_n)

autoplot(model_liniowy_n) + theme_classic()

4.2. Ocena modelu

#dopasowanie modelu
smr_model_lm_n <-summary(model_liniowy_n)

c( R2 = smr_model_lm_n$r.squared, R2_adj = smr_model_lm_n$adj.r.squared)
##        R2    R2_adj 
## 0.8914342 0.8903486
#~89% zmienności zmiennej TOT jest wyjaśniona przez zmienność zmiennej N_BUDYNKI.

# Określenie dowolnego przedziału ufności dla współczynników modelu
confint(model_liniowy_n, level = 0.99)
##                 0.5 %   99.5 %
## (Intercept) -2.413594 3.198192
## N_BUDYNKI    3.198037 3.843294
outlierTest(model_liniowy_n)
##     rstudent unadjusted p-value Bonferroni p
## 65 -6.324628         7.3790e-09   7.5266e-07
## 33  4.742867         7.0854e-06   7.2271e-04
ggplot(pop1km_n, aes(x=TOT, y=N_BUDYNKI)) + 
  geom_point() +
  geom_point(data = pop1km_n[c(65, 33),], aes(x=TOT, y=N_BUDYNKI), color = "red", size = 4) +
  coord_fixed(ratio = 1) + 
  xlab("Liczba ludności") + 
  ylab("Liczba mieszkań") + 
  geom_smooth() + 
  geom_abline(color = 'grey') + 
  theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

#WIZUALIZACJA
par(mfrow = c(2, 4))
plot(model_liniowy)
plot(model_liniowy_n)

meanerror = mean(model_liniowy_n$residuals)
rmse = sqrt(mean((model_liniowy_n$residuals)^2))

print(paste("średni błąd estymacji modelu wynosi", meanerror, "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))
## [1] "średni błąd estymacji modelu wynosi 2.30805383656037e-16 natomiast pierwiastek średniego błędu kwadratowego wynosi 9.14"

5. Wizualizacja wyników modelu na mapie

5.1. Mapa pokazująca oryginalne wartości liczby ludności (TOT)

ggplot() +
  geom_sf(data = pop1km, aes(fill = TOT)) +
  scale_fill_viridis_c() 

  labs(fill = "Liczba ludności", title = "Mapa liczby ludności") +
  theme_minimal()
## NULL

5.2. Mapa pokazującą estymowaną liczbę ludności (EST_POP)

# Przewidywanie EST_POP na podstawie modelu liniowego
pop1km$EST_POP <- predict(model_liniowy, newdata = pop1km)

# Tworzenie mapy z estymowaną liczbą ludności (EST_POP)
ggplot() +
  geom_sf(data = pop1km, aes(fill = EST_POP)) +
  scale_fill_viridis_c() +  # Możesz zmienić paletę kolorów wg potrzeb
  labs(fill = "Estymowana liczba ludności", title = "Mapa estymowanej liczby ludności") +
  theme_minimal()

5.3. Mapa reszt (różnic między wartością obserwowaną TOT oraz EST_POP)

# Obliczenie reszt
pop1km$Reszty <- pop1km$TOT - pop1km$EST_POP

# Tworzenie mapy reszt
ggplot() +
  geom_sf(data = pop1km, aes(fill = Reszty)) +
  scale_fill_viridis_c() +  # Możesz zmienić paletę kolorów wg potrzeb
  labs(fill = "Reszty", title = "Mapa reszt (TOT - EST_POP)") +
  theme_minimal()

5.4. Mapa reszt przeklasyfikowaną do 3 klas (błędy ujemne, brak błędu, błędy dodatnie)

# Przeklasyfikowanie reszt do trzech klas
pop1km$Klasy_reszt <- cut(pop1km$Reszty, breaks = c(-Inf, 0, Inf), labels = c("Przeszacowane", "Niedoszacowane"))

# Utworzenie etykiety dla braku błędu
pop1km$Klasy_reszt[is.na(pop1km$Klasy_reszt)] <- "Brak błędu"
## Warning in `[<-.factor`(`*tmp*`, is.na(pop1km$Klasy_reszt), value =
## structure(c(1L, : niepoprawny poziom czynnika, wygenerowano wartość NA
# Tworzenie mapy przeklasyfikowanych reszt
ggplot() +
  geom_sf(data = pop1km, aes(fill = Klasy_reszt)) +
  scale_fill_manual(values = c("Przeszacowane" = "red", "Brak błędu" = "grey", "Niedoszacowane" = "blue")) +  # Możesz zmienić kolory wg potrzeb
  labs(fill = "Klasy reszt", title = "Mapa przeklasyfikowanych reszt") +
  theme_minimal()

6. Statystyczna i przestrzenna analiza rokładu błędów modelu regresji na podstawie map oraz wyników modelu

7. Mapa rozmieszczenia ludności w siatce 100m, dane pomocnicze - liczba mieszkań.

# Tworzenie interaktywnej mapy rozmieszczenia ludności w siatce 100m
pop_100m_plot <- ggplot() +
  geom_sf(data = pop_100m, aes(fill = TOT)) +
  scale_fill_viridis_c() +  
  labs(fill = "Liczba ludności", title = "Mapa rozmieszczenia ludności w siatce 100m") +
  theme_minimal()

# Konwersja ggplot na obiekt plotly
pop_100m_plotly <- ggplotly(pop_100m_plot)

# Wyświetlenie interaktywnej mapy
pop_100m_plotly